library(SemiParBIVProbit)
## a simulated data example with Bernoulli outcomes
## 
set.seed(0)

n.rep=1  #is the number of replicates

n=1000 #sample size

#d.unob indicates the distribution for the unobsevable, it can #be normal (1),
# t with 4 dof (2), chi2 with 1 dof (3), uniform (-3,3) (4)



sim.func <- function(n.rep=100, n=500, d.unob = 1, beta = 0){



mat.coef <- matrix(NA, n.rep, 9)
true <- rep(NA, n.rep)


for(i in 1:n.rep){

## measured potential confounders
z1 <- runif(n)

## unmeasured confounder
if(d.unob==1){

z2 <- rnorm(n)       
                     
}                     
                     
if(d.unob==2){

z2 <- rt(n, df=4)       
} 

if(d.unob==3){

z2 <- rchisq(n, df=1)       
} 

if(d.unob==4){

z2 <- runif(n, -3,3)       
} 

  

z4 <- runif(n)                     
z3 <- rbinom(n,1,0.5) # IV


## non-linearities
f1   <- function(x) cos(pi*2*x) + sin(pi*x)
f2   <- function(x) x + exp(-30*(x-0.5)^2) 

## treatment assignment

prob.treated <-	plogis(-0.5 + f1(z1) + beta*z2 + 3*z3 + 1.3*z4) 
x <- rbinom(n, 1, prob.treated)

## potential outcomes
p0 <- plogis(-3.5 + f2(z1) + 2*z2 -0.8*z4)                    
p1 <- plogis( 0.5 + f2(z1) + 2*z2 -0.8*z4)                   


y0 <- rbinom(n, 1, p0)
y1 <- rbinom(n, 1, p1)

## observed outcomes
y <- y0
y[x==1] <- y1[x==1] 	     	

table(y)
table(x)
true[i] <- mean(y1)-mean(y0)

## put everything in a data frame
examp.data <- data.frame(z1, z2, z3, z4, x, y)


tr.eq <- x ~ s(z1) + z3 + s(z4)
out.eq <- y ~ x + s(z1) + s(z4)



resProbitProbit <- try(SemiParBIVProbit(list(tr.eq, out.eq), data=examp.data, margins = c("probit","probit")))
if ('try-error' %in% class(resProbitProbit)) next

resLogitLogit <- try(SemiParBIVProbit(list(tr.eq, out.eq), data=examp.data, margins = c("logit","logit") ))
if ('try-error' %in% class(resLogitLogit)) next

resCloglog <- try(SemiParBIVProbit(list(tr.eq, out.eq), data=examp.data, margins = c("cloglog","cloglog")))
if ('try-error' %in% class(resCloglog)) next

resProbitProbitF <- try(SemiParBIVProbit(list(tr.eq, out.eq), data=examp.data, margins = c("probit","probit"), BivD="F" ))
if ('try-error' %in% class(resProbitProbitF)) next


resLogitLogitF <- try(SemiParBIVProbit(list(tr.eq, out.eq), data=examp.data, margins = c("logit","logit"), BivD="F" ))
if ('try-error' %in% class(resLogitLogitF)) next


resCloglogF <- try(SemiParBIVProbit(list(tr.eq, out.eq), data=examp.data, margins = c("cloglog","cloglog"), BivD="F" ))
if ('try-error' %in% class(resCloglogF)) next


mat.coef[i,] <- c( AT(resProbitProbit,   nm.end = "x",type="univariate")$res[2],
                   AT(resProbitProbit,   nm.end = "x",type="simultaneous")$res[2],

                   AT(resLogitLogit,   nm.end = "x",type="univariate")$res[2],
                   AT(resLogitLogit,   nm.end = "x",type="simultaneous")$res[2],

                   AT(resCloglog,    nm.end = "x",type="univariate")$res[2],
                   AT(resCloglog,    nm.end = "x",type="simultaneous")$res[2],

                 
                   AT(resProbitProbitF,  nm.end = "x",type="simultaneous")$res[2], 
 

                   AT(resLogitLogitF,    nm.end = "x",type="simultaneous")$res[2],


                   AT(resCloglogF,  nm.end = "x",type="simultaneous")$res[2]

 

		                   
                     
                     )
                     
                     
                     
print(i)

}

   

true.v <- mean(true)
bias <- abs((colMeans(mat.coef[,1:9], na.rm = T)-true.v)/true.v)*100
names(bias) <- c(         "UnivariateP",
                          "NormalP",

                          "UnivariateL",
                          "NormalL",

		               "UnivariateC",
                          "NormalC", 

		
                          "FrankP",
                          "FrankL",


                          "FrankC" 
   


			  )
			  
			  
			  
mse <- apply((mat.coef[,1:9] - true.v)^2, 2, mean, na.rm=TRUE)
names(mse) <- c  ("UnivariateP",
                          "NormalP",

                          "UnivariateL",
                          "NormalL",

		               "UnivariateC",
                          "NormalC", 

		
                          "FrankP",
                          "FrankL",


                          "FrankC"


			  )


results <- list(bias=bias, mse= mse, mat.coef=mat.coef, true.v= true.v)

results

}

results = sim.func(n.rep=250, n=1000, d.unob = 3, beta = -1.5)

save.image("resminus1.5un3.RData")
bias <-results$bias
round(bias)

rmse <-sqrt(results$mse)
round(rmse,2)







